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Abstract 

Context. Mathematical optimization can be used as a computational tool to obtain the optimal solution to a given problem in a sys- 
tematic and efficient way. For example, in twice-differentiable functions and problems with no constraints, the optimization consists 
of finding the points where the gradient of the objective function is zero and using the Hessian matrix to classify the type of each 
point. Sometimes, however it is impossible to compute these derivatives and other type of techniques must be employed such as the 
steepest descent/ascent method and more sophisticated methods such as those based on the evolutionary algorithms. 
Aims. We present a simple algorithm based on the idea of genetic algorithms (GA) for optimization. We refer to this algorithm as 
AGA (Asexual Genetic Algorithm) and apply it to two kinds of problems: the maximization of a function where classical methods fail 
and model fitting in astronomy. For the latter case, we minimize the chi-square function to estimate the parameters in two examples: 
the orbits of exoplanets by taking a set of radial velocity data, and the spectral energy distribution (SED) observed towards a YSO 
(Young Stellar- Object). 

Methods. The algorithm AGA may also be called genetic, although it differs from standard genetic algorithms in two main aspects: 
a) the initial population is not encoded, and b) the new generations are constructed by asexual reproduction. 

Results. Applying our algorithm in optimizing some complicated functions, we find the global maxima within a few iterations. For 
model fitting to the orbits of exoplanets and the SED of a YSO, we estimate the parameters and their associated errors. 

Key words. Methods:numerical-Stars-55 Cancri-Planets and Satellites:general-ISM:individual(L1448) 



1. Introduction 

' Mathematical optimization can be used as a computational tool 
in deriving the optimal solution for a given problem in a sys- 
tematic and efficient way. The need to search for parameters that 

I cause a function to be extremal occurs in many kinds of opti- 

, mization. The optimization techniques fall in two groups: deter- 
ministic (Horst and Tuy, 119901 ) and stochastic (Guus, Boender 
and Romeijn J1995I I. In the first group, we have the classi- 
cal methods that are useful in finding the optimum solution 

' or unconstrained maxima or minima of continuous and twice- 
differentiable functions. In this case, the optimization consists 
of identifying points where the gradient of the objective func- 
tion is zero and using the Hessian matrix to classify the type of 
each point. For instance, if the Hessian matrix is positive defi- 
nite, the point is a local minimum, if it is negative, the point is 
a local maximum, and if if indefinite, the point is some kind of 
saddle point. However, the classical methods have limited scope 
in practical applications since some involve objective functions 
that are not continuous and/or not differentiable. For these rea- 
sons, it is necessary to develop more advanced techniques that 
belong to the second group. Stochastic models rely on proba- 
bilistic approaches and have only weak theoretical guarantees 
of convergence to the global solution. Some of the most use- 
ful stochastic optimization techniques include: adaptive random 
search (Brooks, 119581 ), clustering methods (Torn, ,19731 ), evo- 
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lutionary computation that includes genetic algorithms, evolu- 
tionary strategies and evolutionary programming (Fogel, Owens 
and Walsh, 1966; Schwefel, [19951 Goldberg, 1989, McCall, 
2005), simulated and quantum annealing (Kirkpatrick, Gellatt 
and Vecchi, il983) , and neural networks (Bounds, [1987]). 

We present a simple algorithm for optimization (finding the 
values of the variables that maximize a function) and model fit- 
ting (finding the values of the model parameters that fit a set 
of data most closely). The algorithm may be called genetic, 
although it differs from standard genetic algorithms (Holland, 
119751 ) in the way that new generations are constructed. Standard 
genetic algorithms involve sexual reproduction, that is, the re- 
production by the union of male and female reproductive in- 
dividuals. Instead, our algorithm uses asexual reproduction, in 
which offspring are produced by a single parent (as in the fission 
of bacterial cells). 

The paper is organized as follows. In Sect. 2, we present and 
describe the main characteristics of our algorithm called AGA 
(Asexual Genetic Algorithm). In Sect. 3, we apply the algo- 
rithm to two kinds of problems: maximization of complicated 
mathematical functions and a model fitting procedure. In the lat- 
ter group, we consider two examples taken from astronomy: a) 
the orbital fitting of exoplanets, and b) the model fitting of the 
Spectral Energy Distribution (SED) observed in a Young Stellar 
Object (YSO). In both cases, we minimize their corresponding 
chi-square function. In Sect. 4, we summarize and discuss the 
results for each case. 
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2. Description of thie AGA (Asexuai Genetic 
Algorithim) 

We consider the problem of finding the absolute maxima of a 
real function of variables i.e., identifying the values of the 
variables (the coordinates of a point in the space of dimen- 
sions) for which the function attains its maximum value. It is 
assumed that the absolute maximum is inside a bounded region 
y where the function is defined. 

Our algorithm proceeds in the following way (see also Fig. 

m-. 

1. Construct a random initial population. The initial popula- 
tion is a set of Nq randomly generated points (in the context 
of evolutionary algorithms, they are also called individuals) 
within the region V. 

2. Calculate the fitness of each individual in the population. The 
fitness is calculated by evaluating the function at each point. 

3. Select a subset of individuals with the highest fitness. Rank 
the points according to the value of the function, and choose 
a subset of A^i points with the highest values of the function. 

4. Construct a new population using the individuals in the sub- 
set. Generate N2 random points within a previously selected 
vicinity E around each of the selected points. 

5. Replace the source population with the new population. The 
new population is the set of (A^i x A^2) + ^1 points that re- 
sults from step 4 plus a clone of each parent. We may choose 
(as we did) A^i and N2 such that (iVi x N2) + Ni = Nq, 
keeping the size of the population Nq unchanged for each 
generation. In this way, one can devise an iterative proce- 
dure. 

6. If the stopping criteria (accuracy, maximum number of gen- 
erations, etc.) have not been met, return to step 2. 

It is clear that this presented algorithm in many aspects resem- 
bles standard genetic algorithms. The key difference is the way 
the new population is constructed: in this version, we propose 
an asexual reproduction with mutation in the sense that the off- 
spring of a parent is a point randomly selected within a nar- 
row neighborhood of the (single) parent point. Hence, the name 
for our algorithm is AGA, meaning Asexual Genetic Algorithm 
(Fig. [Til. We note that in the new generation we always include 
a clone of the parent, that is, we always keep the original seed 
points to be used for the next generation whenever the parents 
are more well suited than their offspring. 

For the vicinity E of each selected point, we used rectangular 
(hyper)boxes of decreasing size. The box around each point was 
centered on the point and had a side length 2Ay along direction 
i in the generation j. In particular, we take, 

A,,^XqP>, (1) 

where 2Aio is the initial length of the box along direction /, and 
/? is a fixed numerical value less than unity (which can be called 
the "convergency factor"). In this way, the length of each box 
side decreases by a factor p in each generation. For instance, if 
we want the side length of the box to decrease by a factor 2 after 
10 generations then, 

P=Qy°= 0.9330 (2) 

Decreasing the size of the vicinity E each generation is in- 
tended to achieve the highest possible accuracy for the position 
of the point at which the function attains its absolute maxima. 



The speed with which the AGA finds a solution depends, of 
course, on the factor p; the lower the value of p, the faster the so- 
lution found. However, if p is too low, the sampling area around 
the points may decrease so fast that the AGA has no time to mi- 
grate to the true solution. On the other hand, if p is too high, 
many of the offspring will never reach the solution within the 
convergency criterion. As a consequence of this, the error in the 
solution will be high and in some cases, the solution will not 
be reached. The adequate value for p depends on the problem 
itself. The optimization of p can be achieved by trial and error. 
However, we have found that a value between 0.4 and 0.6 is ad- 
equate for all the tested problems. 

An alternative way of choosing the box size consists of em- 
ploying the standard deviation of the points contained in the sub- 
set iVi along each dimension. In such a case, the length of the 
sides of the box naturally decreases as the algorithm converges. 
This method works quite efficiently for problems with a few di- 
mensions. Interestingly, the length of the box size in this case 
decreases following a power law such as that in Eq. [1] with a 
moderately high values of p (between 0.5-0.8) for the first few 
generations, abruptly changing to a much lower values (0.2-0.3) 
for the rest of the generations. 

In the case of problems with the large number of parame- 
ters to be estimated, we added an iterative method to the scheme 
presented above. This iterative method consisted of performing 
a series of runs (each one following the scheme shown in Fig.[T]) 
in such a way that the resultant parameters of a run are taken as 
the initial "guess" for the next run. This procedure may be equiv- 
alent to performing a single run allowing for additional genera- 
tions, but this is not the case. The key difference is that for each 
run in the iterative procedure, the sizes of the sampling boxes 
are reset to their initial values, i. e., each run starts searching for 
solution using boxes of the same size as those used in the first 
run but centered on improved initial values. We consider to have 
reached the optimal solution when the values of the parameters 
do not change considerably (within a tolerance limit) after sev- 
eral iterations, and the found to have reached a limiting 
value (see sect. 3.2.1 and Fig.|2]for an example). 

We find that this iterative strategy guarantees the conver- 
gence to the optimal solution, since it avoids the potential danger 
of using values of p that do not allow the AGA to drift (mi- 
grate) to the "true" solution within a single run. Problems that 
involve the finding of a large number of parameters are poten- 
tially subject to this risk, since each parameter may require a 
different value of p. Furthermore, problems with a large number 
of parameters become particularly difficult when the values of 
the different parameters differ by several orders of magnitude, 
as in the fitting of the orbits of exoplanets and the fitting of the 
SED of YSOs (see Examples 1 and 2 in Sect. 3). In this case, the 
iterative method has proven to be particularly useful; the solution 
usually improves considerably after several iterations. 

In the following section, we describe some applications of 
AGA. We have divided the applications into two groups depend- 
ing on the type of problem to solve; the maximization of com- 
plicated functions and model fitting in astronomy. 

3. Applications of thie AGA 

We separate the optimization problems in two groups. In the 
first group, we consider functions of two variables where many 
classical optimization methods present formidable difficulties 
in finding the global maximum. In the second group, we show 
two typical examples of model fitting in astronomy that can be 
treated as minimization procedures. 
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Figure 1. Basic diagram for the implementation of the Asexual 
Genetic Algorithm (AG A). First we generate a random initial 
population. Then, we evaluate the fitness of each individual in 
this population and select those which have the highest fitness. 
A new generation is constructed by an asexual reproduction (see 
text) which replaces the older one. If the stopping criteria are 
met, we stop. If not, we use these individuals as an initial popu- 
lation and start again. 



3. 1. Optimization of functions of two variables 

There are many examples of functions that are not easy to op- 
timize with classical techniques such as the simplex method, 
the gradient (or Newton-Raphson) method, the steepest ascent 
method (Everitt, I1987I I. among others. In such cases, the exis- 
tence of many maxima (or minima) and the sharpness of the 
peaks can represent a serious problem. Because of this, the 
standard Genetic Algorithms or GA are successfully applied in 
searching for the optimal solution. We consider some typical ex- 
amples treated by this technique, which are shown below. 

3.1 .1 . Example 1 

We consider the following function (Charbonneau, ! 1 995T l : 

/ {x, y) — [IQx {1 — x)y [1 ~ y) sin {mrx) sin {mry)]^ (3) 

where the variables x, y G [0, 1] and n G N — {0}. Identifying 
the global maximum of this function for large n is a difficult task 
because there are many local maxima that differ little in value 
but are separated by deep "valleys" in the two-dimensional land- 
scape. Techniques such as the steepest ascent/descent and the 
conjugate gradient method are local methods that work well if 
/ {x, y) is a smooth function that can be differentiated at least 
once and a single maximum exists in the domain under consid- 
eration. 



Table 1. Values of the initial parameters used in AGA for the 
four examples presented in this work. 



— -____AGA's parameter 
Example — — 


P 




N2 


iVo 


Bi-dimensional functions (2) 


0.6 


10 
20 


9 
19 


100 
400 


Extrasolar planets 
YSO 


0.45 


30 


29 


900 



In Fig. |2] we show the solution to this problem (n = 9) 
by applying AGA after fixing the initial population size, num- 
ber of parents, number of descendants, and convergency factor 
(Table [1]). Our graphs have the same format as those presented 
by Charbonneau ( |1995l l to facilitate their comparison. 

In the first panel, we start with a population of A'o = 100 
individuals, i. e., a set formed by 100 random points represent- 
ing the candidate solutions to the global maximum distributed 
more or less uniformly in parameter space. After 5 generations 
(second panel), a clustering at the second, third and fourth max- 
ima is clearly apparent. After 10 generations, the solutions al- 
ready cluster around the main maxima. At the 25th generation, 
we have reached the maximum in {x*,y*) — (0.5,0.5) with 
/ {x*,y*) = 1 with an accuracy of ^ lO"*". We note that at the 
25th generation all the 100 individuals have reached the maxi- 
mum with at least this accuracy. 

In the last panel, we show the evolution in the fittest "phe- 
notype " with the number of generations as plotted for two sizes 
of the population, A^o = 100 and A^o = 400; in other words, 
we measure the deviation in the function value for the maxima 
points identified in each generation from the "true" maximum, 
that is, / (0.5*; 0.5*) — 1. It is evident from Fig.|2]that a larger 
size of the population causes the maximum to be reached in a 
lower number of generations. 

In Fig. |3] we show the evolution in the global solution using 
AGA and the results obtained by Charbonneau (1995) using a 
traditional GA. In the first panel for both works, an initial popu- 
lation of 100 individuals is randomly selected. After 10 genera- 
tions, all the individuals in our algorithm have clustered around 
the maxima; in contrast, the GA has started to form some clus- 
ters. In the third panel, after 20 generations, AGA has found the 
global maxima, while the GA continues to search among clus- 
ters for a global solution. In the fourth panel, the convergence 
to the global solution is shown for both works, our AGA hav- 
ing reached the solution in 25 generations with higher accuracy 
than the 90 generations employed by the GA. In the fifth panel, 
we finally can compare the evolution in the fittest phenotype for 
each generation, calculated to be 1 — / {xmax,ymax)- For AGA, 
these differences are much smaller than those obtained with the 
GA method, which indicates that our algorithm is more accu- 
rate than GA. However, we note that the solution presented in 
Charbonneau's work was limited by the number of digits used 
in the encoding scheme. The GA algorithm can also reach high 
levels of precision by changing the encoding, but at the cost of 
slower convergence, which would make the GA algorithm be- 
come even more computationally expensive than the AGA. The 
lower number of generations required by AGA to reach a global 
solution demonstrates that AGA is more efficient than the GA 
method. 



4 



J. Canto, S. Curiel and E. Martmez-Gomez: A simple algorithm for optimization and model fitting 



0.5 



N = 

T' ' i » .. '• 

"o (^p«(^|©@ o . _ 

mmw V 

: Oj@©@l©(S) o •; 



N = 5 



O (o) ( 



N = 10 

) (o) O 



0.5 



-o (( 



N = 20 











© ® 




_ 












- ■ c 






_ 




- 




O® © © o - 




■ T 






Figure 2. A solution to the model optimization problem based on the idea of the Asexual Genetic Algorithm. The first five panels 
show the elevation contours of constant / (Eq. [3] with n = 9) and the population distribution of candidate solutions (each one 
contains 100 points), starting with the initial random population (in the standard GA it is defined as the "genotype") in (a) and 
proceeding on through the 25th generation on (e). In the sixth panel, we show the evolution of the fittest "phenotype " assuming two 
sizes of the population, A^o = 100 and A^o = 400. 



3.1.2. Example 2 

The following example is proposed as a function test in the 
PIKAIA's user guide (Charbonneau and Knapp, I1995I ). The 
problem consists of locating the global maximum of the func- 
tion: 



/ (r) = / (x, y) = cos'^ (tnrr) exp 



(4) 



[1 
2a^ J 

where a and b are known constants, r is the radial distance 

given by the expression r = sj {x — xq)^ + {u ^ Uo)'^^ c is the 
width of the Gaussian and the position of the maximum peak is 
{xq, yo). Observed from above, this function appears like con- 
centric rings of similar widths and amplitudes (see FigHJi. The 
difficulty in finding the maximum peak of this function is that 
the area of the concentric rings is much larger than the area at the 
center containing the maximum peak (see FigHJi. This property 
ensures that most of the algorithms searching for the maximum 
peak fail because they usually become fixed inside in one of the 
rings and identify a local maximum instead. We experimented 
with different values of a, b, and (a = 2, 4; 6 = 3, 9 and 
tr^ = 1, 2, 4), as well as different centers (xq, yo). In all cases, 
using AGA with the parameters shown in Table [H the "true" 
maximum peak was found in just a few steps, typically in less 
than 100 generations, and with an accuracy of 10^*^. The results 



Table 2. Values of the parameters a, b, and assumed for 
the function / (r) = cos"- (birr) exp (^—j^^ and the number of 
generations needed to reach an error tolerance of 10^^. 



a 


b 


a' 


Ngen 


2 


3 


1 


65 


2 


3 


2 


69 


2 


3 


4 


66 


2 


9 


1 


94 


2 


9 


2 


94 


2 


9 


4 


91 


4 


3 


1 


69 


4 


3 


2 


66 


4 


3 


4 


58 


4 


9 


1 


94 


4 


9 


2 


101 


4 


9 


4 


98 



are summarized in Table |2] Changing the position of the center 
did not substantially change the number of generations needed 
to reach the selected tolerance level. 
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Figure 3. Comparison of our results with those obtained by Charbonneau (TWS"). We have maintained the same format and removed 
the curve of the evolution of the median-fitness individual in the fifth panel of Charbonneau's work to facilitate direct comparison. 
In the first panel we start with an initial population of 100 individuals and both algorithms start to search the global maxima. While 
the GA employs 90 generations to reach the solution, our AGA just requires 25 (fourth panels). The fifth panels show the evolution 
of the fittest phenotype with the number of generations for AGA (left) and GA (right). Note that indeed, AGA reaches an accuracy 
of ^ 10^^ in less that 30 generations while GA attains an accuracy of ^ 10^^ in 100 generations. 



3.2. The parameter estimation problem in astronomical 
models 

In the physical sciences, curve or model fitting is essentially an 
optimization problem. Giving a discrete set of N data points 
{xi,yi) with associated measurement errors a-i, one seeks the 
best possible model (in other words, the closest fit) for these data 
using a specific form of the fitting function, y (x). This function 
has, in general, several adjustable parameters, whose values are 
obtained by minimizing a "merit function", which measures the 
agreement between the data and the model function y (x). 

We suppose that each data point yi has a measurement 
error that is independently random and distributed as a nor- 
mal distribution about the "true" model with standard deviation 
(Ti . The maximum likelihood estimate of the model parameters 
(ci, Cfc) is then obtained by minimizing the function, 

2 _ V" - y (^»; ci, ■■,cfe)y 
^ = ^ [ J 

3.2.1 . Evaluation of the error estimation 

The experimental or observational data are subject to measure- 
ment error, thus it is desirable to estimate the parameters in the 
chosen model and their errors. In the straight-line-data fitting 
and the general linear least squares,we can compute the stan- 
dard deviations or variances of individual parameters through 
simple analytic formulae (Press et al. . 119961 1. However, when we 



attempt to minimize a function such as Eq.|5] we have no expres- 
sion for calculating the error in each parameter. A good approach 
to solve this problem consists of building "synthetic data sets". 
The procedure is to draw random numbers from appropriate dis- 
tributions so as to mimic our clearest understanding of the un- 
derlying process and measurement errors in our apparatus. With 
these random selections, we compile data sets with exactly the 
same numbers of measured points and precisely the same val- 
ues of all control or independent variables, as our true data set. 
In other words, when the experiment or observation cannot be 
repeated, we simulate the results we have obtained. 

We compiled the synthetic data with the process illustrated 
in Fig. |5]and from the following expression: 



y, ^y, + a, (2^ - 1) (6) 

where y^ represents the i-th data point of the new set, y i is the i-th 
data point of the original set, ai is the error associated with the i- 
th data point, and ^ is a random number within the interval [0,1]. 
Using Eq. |6] we generated synthetic data sets with the original 
errors because they are related to the measurements. 

For each of these data sets, we found a corresponding set of pa- 
rameters (see Fig. |5]). Finally, we calculated the average values 
of each parameter and the corresponding standard deviation, as 
the estimates of the parameters and their associated errors, re- 
spectively. We applied this procedure to the following examples. 
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Figure 4. Profile of the bi-dimensional positive cosine function. In the upper panels, we take a = 2, b = 3 and = 2; in the 
lower panels, we take a = 4, 6 = 9 and cr^ = 4. In both cases, the position of the maximum peak is {x, y) — (0.25, 0.25), where 
/(0.25, 0.25) = 1. 
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Figure 5. From the original set of measurements/observations, 
several synthetic data sets are constructed. This is achieved by 
adding to the dependant variable a random number whose abso- 
lute value is within the estimated errors of the original data. For 
each of these synthetic data sets, a set of fitted parameters is ob- 
tained. The average value of each parameter and the correspond- 
ing standard deviation are taken as estimates of the parameters 
and their associated errors. 



3.2.2. Fitting tine orbits of Extrasolar Giant Planets 

We use the merit function given in Eq. |5] and the algorithm de- 
scribed in Sect. 2 to solve an interesting and challenging task in 
astronomy: curve or model fitting of data for the orbits of the 
extrasolar planets. 



The first extrasolar planet was discovered in 1995 by Mayor 
and Queloz (1995) and, according to The Extrasolar Planets 
Encyclopaedi^l until February 2009 there had been 342 can- 
didates detected. Most of them (316) were revealed by radial 
velocity or astrometry of 269 host stars, that is, by Keplerian 
Doppler shifts in their host stars. Doppler detectability favors 
high masses and small orbits depending mainly on the present 
Doppler errors achievable with available instruments. In addi- 
tion, the precision of the Doppler technique is probably about 
3 m s owing to the intrinsic stability limit of stellar photo- 
spheres. This technique is sensitive to companions that induce 
reflex stellar velocities, K > 10ms~^, and exhibit orbital pe- 
riods ranging from a few days to several years, the maximum 
detectable orbital period being set by the time baseline of the 
Doppler observations. The remaining exoplanets were detected 
by other techniques: microlensing (8), imaging (11), and timing 
(7). For this example we only refer to the planets detected by 
radial velocity. 

We now consider a system consisting of a central star of mass 
Af,, suiTounded by Np planets in bounded orbits. Assuming that 
the orbits are unperturbed Kepler orbits, the line-of-sight veloc- 
ity of the star relative to the observer (Lee and Peale, |2003T l is, 

v^, = vq + ^ Kj [cos {fj +ujj) + ejcosujj] (7) 



http://exoplanet.eu/catalog 
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where ej is the eccentricity, ojj is the argument of the pericenter, 
fj is the true anomaly, Kj is the velocity amplitude of planet j, 
and vq is the line-of-sight velocity of the center-of-mass relative 
to the observer. The true anomaly is related to time by means of 
Kepler's equation. 



E 



■sin (E) = ^{t-toj) 



and 



tan 



1 



-tan 



(8) 



(9) 



where Tj and toj are the period and time, respectively, of peri- 
center passage for planet / Thus, for each planet, there are 5 free 
parameters: ej , Tj, toj, Kj and cuj. Additionally, there is the 
systemic velocity vq. In total, we have (5 x Np) + 1 free param- 
eters, which are simultaneously fitted using AGA. 

When these basic parameters for each planet are known 
through model fitting, we can estimate the semi-major axes of 
the orbits, Uj, and the masses of each planet, mj. To do this, we 
have to make other simplifying assumptions. For instance, in the 
simplest model of totally independent planets. 



nijsin (ij) 



1 



(10) 



from which we can determine either nij by assuming sin ij — 1 
or nijsinij by neglecting mj in front of the mass of the star, 
M*. We can then estimate aj using Kepler's third law. 



rp2 



G {Ah + rrij) 



(11) 



where G denotes the gravitational constant. 

We applied the AGA presented in Sect. 2 to the measured 
radial velocities of the main-sequence star 55 Cancri, published 
by Fischer et al. (|2008] l. 

The data set contains 250 measurements completed at the 
Lick Observatory from 1989 to 2007, and 70 measurements 
made at the Keck Observatory from 2002 to 2007 of the ve- 
locity of 55 Cancri. Data from Lick has measurement errors 
of ~ lOms-i (1989-1994) and 3-5 ms"! (1995-2007). Data 
from Keck has measurement errors of 3 m s^^ for data acquired 
prior to 2004 August and 1.0-1.5 ms^^ thereafter. The mea- 
sured data have a maximum amplitude of ^ 150 ms^^ and are 
of excellent signal-to-noise ratio. The orbital parameters were 
established by the detection of four planets in all Doppler mea- 
surements. To test a possible stellar systemic velocity remnant 
in the data, we added an additional parameter to the systemic 
velocity of the star obtaining a value of 17.28 ms~^. We fit- 
ted the orbits of 1, 2, 3, and 4 planets to this data (see Fig. 
|6l). This exercise shows that the fitting solution improves when 
considering the orbit of more planets. The rms of the residuals 

and the Xred, [xled = jVp„.„t,-jVpl„„et,„-i ) values, improve 
from values of 33.5 and 9.72, respectively, when fitting the orbit 
of one planet, to values of 7.99 and 2.03, respectively, when fit- 
ting the orbit of four planets (see TableO. In general, the values 
of the rms, VX^ed ^'-^ ™^ '-^^ derived parameters compare well 
with those obtained by Fischer et. al. ( I2008I I. In Tables [3] and |4] 
we summarize our results. 



Table 3. Our values of the rms of the residuals and the fits 
for the orbital fitting problem using AGA. For comparison, in the 
last two columns we show their corresponding values obtained 
by Fischer et al. ( |2008] l. 



1 

Planet 


This work 
rms (m s ') 


This work 


Fischer et al. 
rms (m s ') 


Fischer et al. 


b 


33.5 


9.72 


39 


10 


b, c 


10.69 


3.38 


11.28 


3.42 


b, c, d 


9.69 


2.49 


8.62 


2.50 


b, c, d, e 


7.99 


2.03 


7.87 


2.12 



'The single-planet fitting was obtained in Marcy et al. ( |2002i l. 



With the exceptions of the third and fourth planet fittings, 
our rms values are lower than those obtained by Fischer et al. 
(12008'). Our values are however lower in all cases. 

From the parameters shown in Table|4]we were able to derive 
the mass of the planet, Mpsini (Mj), and the major semiaxis a 
(AU). For the four-planet fitting, we found (Table|5]l that our val- 
ues compare well with those reported by Fischer et al. ( I2008I I in 
their five-planet model. Except for the first planet, all the values 
for the major semiaxis have a standard eiTor lower than those 
obtained in the Fischer's model. In the case of the mass of each 
planet, our values have smaller standard errors. 

We included the iterative scheme discussed in sect. 2 in 
fitting the orbit of the planets. In Fig. |7] we show the value of 
the x^g^ for the four-planet fit calculated after each run. In this 
case, we started with x^ed — 4.16 and after the 10th iterative 
run, the x^ed diminished to 4.15, which means that we had 
not found the optimal solution. During the first iterative runs, 
we note that the value of x^ed decreased rapidly. After about 
100 iterative runs, the Xred ^'^^ changed significally and 
converged to a fixed value (there are only slight variations in the 
last decimals). At this point, we can be assured that the value of 
the parameters have been found. Finally, we estimated the errors 
in fitting the orbits of the four planets as described in sect. 3.2.1. 
We present the estimated errors in Tables |4] and |5] 



3.2.3. Fitting tine spectral energy distribution (SED) for a 
YSO 

A spectral energy distribution (SED) is a plot of flux, bright- 
ness or flux density versus frequency/wavelength of light. It is 
widely used to characterize astronomical sources i. e., to iden- 
tify the type of source (star, galaxy, circumstellar disk) that pro- 
duces these fluxes or brightness. Modelling the observed SEDs 
can help us to infer the temperature and size, among the other 
physical parameters of the source. As examples in radio astron- 
omy, a SED with a negative spectral index ~ —0.7 would indi- 
cate the presence of a synchrotron radiation source; in infrared 
astronomy, SEDs can be used to classify T-Tauri stars; in galactic 
astronomy, the analysis of the SEDs leads to the determination 
of the respective roles of the old and young stellar populations 
in dust-grain heating. 

For the reasons explained above, it is interesting to find the 
adequate model for an observed SED. We consider the observa- 
tions reported by Curiel et al. ( |2009l l in the L1448 region. This 
cloud is part of the Perseus molecular cloud complex located at 
a distance of ~ 300 pc. We are interested in fitting the observed 
SED for a couple of reasons: it contains an extremely young 
and highly collimated bipolar outflow (Bachiller et al. |19901 l and 
seems to be a site of very recent star formation based on some 




Figure 6. The velocities and fits for each of the four planets are shown separately for clarity by subtracting the effects of the 
other planets. That is, for planet labelled b, planets c, d and e have been removed from the data using the parameters found in the 
simultaneous four-planet fitting. For planet labelled c we have removed planets b, d and e. For planet labelled d we have removed 
planets b, c and e. For planet labelled e we have removed planets b, c and d. We have added the fitted systemic velocity for the star, 
which value is 17.2826 m s"'. The last panel (bottom and right) shows the residuals between the data and the model. The first (up 
and left) panel shows the raw data. 



Table 4. Our estimated values for the five parameters of the four exoplanets around 55 Cancri. The planets are listed in order of 
increasing orbital period, and the planet designations, b-e, correspond to the notation given by Marcy et al. ( I2002I I and Fischer et 
al. (i2008i) . The value for the \/x^ is also included. 



Planet 


T (days) 


to (JD) 


e 


Lu (deg) 


Ki (m s ') 




e 


2.8170 ± 0.0932 


10000.0046 ± 0.3250 


0.07 ± 0.0016 


250.2326 ± 0.4613 


5.4311 ±0.2241 




b 


14.6515 ±0.0002 


10002.8917 ± 0.0549 


0.0145 ± 0.0005 


130.9176 ±0.4591 


71.7606 ±0.3140 


2.0314 


c 


44.3298 ± 0.0059 


9989.9237 ± 0.4637 


0.0853 ± 0.0014 


78.2384 ± 0.4662 


9.9820 ± 0.2075 


d 


5218.3339 ± 0.5246 


12500.7572 ± 0.3713 


0.0250 ± 0.0006 


180.2123 ±0.2645 


46.6872 ±0.1431 





Table 5. Our derived values for the mass of the planet and the major semiaxis for the four-planet fitting. For comparison, we also 
show the values obtained by Fischer et al. (i2008i) for their five-planet model, where we have removed the fifth body. 



This work This work Fischer et al. Fischer et al. 

Mpsini (Mi) a (AU) Mpsini {Mj) a (AU) 

e 0.0361±0.0014 0.0383±0.0008 0.034±0.0036 0.038±1.0 x 10"" 
b 0.8285±0.0036 0.1148±1.0 x 10"'' 0.824±0.007 0.115±1.1 x 10"*^ 
c 0.1661±0.0035 0.2402±2.1 x 10"^ 0.169±0.008 0.240±4.5 x 10"^ 
d 3.8201±0.0117 5.7705±0.0004 3.835±0.08 5.77±0.11 



observations (Anglada et al. 1989; Curiel et al.'1990' Barsony et 
al.[i998; Girart & Acord 2001 ; Reipurth et al. 2002). 

The data are taken from Curiel et al. (|2009l l, who carried 
out a fit of the SED by assuming that there are three main com- 
ponents contributing to the flux at different wavelengths. Thus, 



the fitting function was given by the contribution of an optically 
free-free component and two grey bodies; 
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V 



Figure 7. Convergence of the Xred value with the number of runs 
(do not confuse with the number of generations done within 
AGA). The curve shows how Xred converges to a "limiting" 
value in the sense that the variations occur in the last decimals 
(vertical line). As the number of runs increases the value of Xred 
diminishes which means that we have not reached the optimal 
solution. Only after about 100 runs the Xred value has converged 
and at this point the parameters do not significally change after 
subsecuent runs. 



Si, — ClLu'^ 



C3 1 



1 



(12) 

where w = ^ is the frequency normalized to a reference fre- 
quency, 1^0 ■ 

The free parameters to be estimated, by the minimization of the 
chi-square (Eq. |5]l, are identified as ci, C2, ciq. They are de- 
fined in terms of their corresponding physical parameters as: 
ci = Fc is the flux of optically thin emission with a spectral 
index given by C2 = a, which corresponds to the free-free emis- 
sion coming from a thermal jet; — Fi = HMl^ii. jg ^ refer- 
ence flux; C4 = Ti is the dust opacity evaluated at the refer- 
ence frequency i/q; cs = /3i is the dust emissivity index; and 

C6 



~ is related to the temperature. The last four parame- 



feTi 

ters correspond to the emission originating in the first grey body 
(probably a molecular envelope). Finally, we propose that a sec- 
ond grey body exists that is a circumstellar disk surrounding the 
young protostar. Similarly, the parameters cy, cs, cg, and cio are 
related to the physical parameters of the second grey body. 
In the previous expressions, ly is the frequency at which the 
source is observed, h is the Planck's constant, il is the solid an- 
gle subtended by the source, c is the speed of light, r is the op- 
tical depth at the frequency t^o for each component, f3 is the dust 
emissivity index for each component, and k is the Boltzmann's 
constant. We assumed the characteristic frequency of the source 
to be 1^0 — — (Curiel et. al., ,2009 ). Table |6] summarizes the 
values of the fitted parameters in the model described by Eq. [12] 
of the SED for a YSO in the L1448 region. Table|6]also includes 
the estimated errors in the fitted parameters, as described in sect. 
3.3.1. Figure [8] shows the observed data together with the best 
fit. 



Table 6. Estimated values for the ten parameters in the model 
of the SED for a YSO in the L1448 region. 



Parameter 


Value 




Fc (mJy) 


0.9295 ± 0.0330 




a 


-0.0169 ± 0.0166 




Fi (mJy) 


1.7109± 0.1141 






0.0728± 0.0062 






0.8999 ± 0.0575 


1.9720 


Tdust 1 (K) 


75.7522 ± 0.5273 


F2 (mJy) 


0.0068 ± 0.0008 




T2 


34.9140 ± 40.7209 




P2 


1.5992 ±0.3254 




Tdust 2 (K) 


330.8240 ± 1.2946 






Figures. The Spectral Energy Distribution (SED) for a YSO 
in the L1448 region. The dots represent the observations and 
the bars their associated measurement errors. The straight line 
corresponds to the continuum flux, the dashed line is the first 
grey body, the dotted line is the second grey body and the solid 
line is the sum of the three contributions. 



From Fig. [8] we can conclude that the fitted model with the 
parameters shown in Table |6] is adequate for explaining the ob- 
servations of the YSO in the L1448 region. The fitted curve lies 
within the measurement error bars (except one point located at 
high frequencies) and the SED consists of a continuum compo- 
nent and the emission from two grey bodies. Further astrophys- 
ical implications of these results can be found in Curiel et al. 
( .2009,) . 



4. Summary and Conclusions 

We have presented a simple algorithm based on the idea of ge- 
netic algorithms to optimize functions. Our algorithm differs 
from standard genetic algorithms mainly in mainly two respects: 
1) we do not encode the initial information (that is, the initial set 
of possible solutions to the optimization problem) into a string 
of binary numbers and 2) we propose an asexual reproduction as 
a means of obtaining new "individuals" (or candidate solutions) 
for each generation. 
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We have then applied the algorithm in solving two types of 
optimization problems: 1) finding the global maximum in func- 
tions of two variables, where the typical techniques fail, and 2) 
parameter estimation in astronomy by the minimization of the 
chi-square. For the latter case, we considered two examples: fit- 
ting the orbits of extrasolar planets associated with the star 55 
Cancri and fitting the SED of a YSO for the L1448 region. 

We found that our algorithm has several advantages: 

- It is easy to implement in any computer because it does 
not require an encoding/decoding routine, and the new 
generations are constructed by the asexual reproduction of 
a selected subset (with the highest fitness) of the previous 
population. 

- The algorithm does not require the evaluation of standard 
genetic operations such as crossover and mutation. This is 
replaced by a set of sampling rules, which simplifies the 
creation of new generations and speeds up the finding of the 
best solution. 

- When the initial "guess" is far from the solution, AGA is 
capable to migrate searching for the "true" optimal solution. 
The final solution is usually achieved after a few hundred 
generations, in some cases, even faster. 

- In some difficult cases, such as the fitting of the orbits 
of several exoplanets and the fitting of the SED of YSOs 
(with several components), AGA finds the solution in a few 
hundred generations. In these cases, an iterative scheme 
(several runs of AGA using the solutions of one run as the 
initial "guess" of the next run) can help to improve finding 
the "true" solution. This is particularly useful when the 
variables have different orders of magnitude, which causes 
the different variables to not converge to the final solution at 
the same time. 

- As a consequence of the previous points, AGA becomes 
computationally less expensive than the standard version 
(GA). The convergence of AGA is reached in just a few gen- 
erations. 

From the example of the orbital fitting for exoplanets around 
55 Cancri, we can conclude that our algorithm gives parameters 
that compare well with those obtained by Fischer et al. (I2008I I. 
Marcy et al. (2002) also used the standard GA to obtain the pa- 
rameters for a third planet, and they concluded that the GA fails 
to significally improve the value. 

In the case of the SED for a YSO in the L1448 region, we 
find that the most adequate model contains the contribution of 
two grey bodies at high frequencies and free-free continuum 
emission, through a power law, at low frequencies. 

Generally speaking the implementation of any kind of GA is 
advantageous in the sense that we do not have to compute any 
derivative of the selected fitness function allowing us to optimize 
functions with several local maxima points. 

In this work, we have also implemented a method to estimate 
the error associated with each parameter based on the generation 
of "synthetic data sets". This method is easy to implement in 
any type of measured/observed data and constitutes a way of 
estimating the errors in the parameters when the minimization 
of the is employed. 
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